Investigations of a Two-Phase Fluid Model 



B. T. Nadiga 
Theoretical Division, Los Alamos National Lab., 
Los Alamos, NM 87545 

S. Zaleski 

Laboratoire de Modelisation en Mecanique, URA 229 CNRS 
Univ. P. et M. Curie (Paris 6), 4 place Jussieu 
75252 Paris Cedex 05, France 

(Submitted to European Journal of Mechanics B: Fluids) 



Abstract 

We study an interface-capturing two-phase fluid model in which the interfacial tension is mod- 
elled as a volumetric stress. Since these stresses are obtainable from a Van der Waals-Cahn-Hilliard 
free energy, the model is, to a certain degree, thermodynamically realistic. Thermal fluctuations are 
not considered presently for reasons of simplicity. The utility of the model lies in its momentum- 
conservative representation of surface tension and the simplicity of its numerical implementation 
resulting from the volumetric modelling of the interfacial dynamics. After validation of the model 
in two spatial dimensions, two prototypical applications — instability of an initially high- Reynolds- 
number liquid jet in the gaseous phase and spinodal decomposition in a liquid-gas system — are 
presented. 
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1. Introduction 

The complexity of the physics involved at interfaces precludes any simple but comprehensive 
mathematical model of it and thus the proper simulation of interfaces is a difficult problem both 
from the physical and the computational points of view. Nonetheless, in several circumstances, it 
is useful to have a model of interfaces that allows for a degree of microscopic realism without being 
constrained by the full details that arise when interface physics are modelled in a way consistent 
with all the basic principles of thermodynamics and particle mechanics. Thus the utility of models 
that capture only some important features of the real world, but are computationally easy to deal 
with. 

Our model is built by adding to the Navier-Stokes equations a stress tensor that is derived 
from the van der Waals-Cahn-Hilliard free energy and results in spontaneous phase separation 
and surface tension. There is a single species of particles, of density p, with the usual continuity 
equation. The van der Waals equation of state ensures the presence of two basic states, a "liquid 
state" of density pi and a "gas" state of density pq. For simplicity, the temperature in the 
equation of state is fixed. While this assumption is physically justifiable in specific flow situations, 
it simplifies the description of the flow by immediately decoupling the energy equation from the 
mass and momentum equations. 

An important feature of our model is that there are special steady state solutions of the model 
mass and momentum conservation equations that connect a gas phase to a liquid phase through a 
smoothly varying density profile of thickness £ . Such an interface is an equilibrium state of the 
Cahn-Hilliard model and this ensures that we describe scales smaller than the physical thickness 
of the interfaces with some degree of realism. Because our model tends to create and maintain 
such a phase separated state, two fundamental physical phenomena are captured: (i) When initial 
conditions involve only thin interfaces (i.e., the hydrodynamic scales, radii of curvature etc... are 
much larger than £), and densities are either pi or pq away from interfaces, the model simulates 
a gas-liquid flow with a free interface between the two phases. In this regime the flow approaches 
real incompressible two-phase flow when the Mach number is small, (ii) When initial conditions 
impose densities in the unstable region the phases will spontaneously separate. (It should be noted 
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that in the early stages of spinodal decomposition the phase separating flow is far from being 
incompressible.) 

Our motivation to study this kind of a model lies in (i) the existence of a spinodal decomposition 
regime, (ii) the fact that the model lends itself to a simple and robust numerical implementation, 
and (iii) that it is a stepping stone to the more complete thermodynamic model summarized in 
the Appendix. Interestingly, the first two characteristics of the model are shared in common with 
the immiscible lattice-Boltzmann models [Appert et al. 1995; Alexander, Chen and Grunau, 1993], 
which are themselves not thermodynamically correct. We note that, unlike most lattice-Boltzmann 
models (with the exception of the recent model by Swift, Osborn, and Yeomans, (1995)), the stresses 
inside an equilibrium interface (thus at uniform temperature) are correct in our model. Further, 
in lattice-Boltzmann models, the tangential velocities are not continuous when pi ^ po, [see 
Rothman and Zaleski, 1994, Ginzburg 1994]. Thus, we may say that our model, while preserving 
the advantages of the lattice-Boltzmann models, eliminates some of its drawbacks. 

Another analogy may be drawn with some surface tension schemes used in other interface 
simulation methods. As in [Brackbill et al. 1993] surface tension is made easier to represent 
numerically because the surface tension stresses are spread continuously over an interface region of 
finite thickness. Moreover, as in [Lafaurie et al. 1994] the surface tension terms conserve momentum 
exactly. 

In section 2, we present the model, interpret the momentum flux components specific to the 
interface and obtain the surface tension coefficient. In section 3 we describe briefly the numerical 
method. In section 4, we present some numerical simulations of the model: we show first that 
the equilibrium densities and pressures on the liquid-gas interface are predicted by the model in 
accordance with the Maxwell equal area rule. We also show that Laplace's law is recovered. Next 
we present the simulation of the instability of a high Reynolds number liquid jet in the gaseous 
phase at two different values of the surface tension. Finally an example simulation of isothermal 
spinodal decomposition with the model is presented. 
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2. The model and the Korteweg interfacial stresses 

We start with the mass and momentum conservation equations: 

d t p + dj(puj) = 0, (1) 
dt(pui) + dj(puiUj) = dj(7ij. (2) 

We decompose the stresses <7 8 j as <7 8 j = —po^ij + + where po is the bulk pressure assumed 
to depend only on the density p through an equation of state po = po(p), and is the viscous 
stress tensor 

a\f = p i^djU, + d t Uj - • uS tJ 

The stresses aff depend on the density gradient and is the key feature of our model. These stresses 
may either be postulated in the manner of Korteweg (1901) or derived thermodynamically from the 
Van der Waals-Cahn-Hilliard free energy (see references in the Appendix as well as the historical 
review of D. D. Joseph (1990)) 



i|Vp)| 2 + P V 2 p) k 



d t pdjp 



(4) 



where A is a constant parameter that controls the strength of the surface tension effect. While sev- 
eral equations of state po(p) will serve our purpose, for definiteness, and in the numerical simulation 
that follows, we choose a van der Waals equation 



where 

Pc = ^PcRgOc, p' = —, d ' = TT 
o p c v c 

and R g = R/M is the universal gas constant divided by the molecular weight and the subscript c 
refers to the critical values (see also the Appendix, Eq. (15)). This equation of state allows for the 
existence of two phases of different densities: there is a range of densities over which the modeled 
fluid is mechanically unstable, i.e., an increase in density results in a decrease in pressure, leading 
to a separation into two different phases. Unlike gas-dynamical models in which the mass and 
momentum balance equations must be supplemented with an equation for either energy or entropy, 
an "artificially compressible" model such as ours needs no further equation. Our model is thus 
entirely described by Eqs. (l)-(5). 



5- 



Considering the stress tensor 




(6) 



where 



T = 



( 



(d x 2 P - d y 2 p)/2 

dxpdyp 



d x pdyp 

(d X 2 P - dy 2 



) 



The first term in Eq. (6) can be thought of as a pressure term acting to smooth gradients in the 
density field, with no contributions to the surface tension. The second term Tj is traceless and 
acts to simulate the surface tension: consider a planar interface along the a; — direction. The only 
nonzero gradients are along the y— direction and we find that 



The mechanical force per unit length along the y-direction is A {T yy — T xx )dy. Identifying it 
with the surface tension coefficient a gives us 



where a is the surface force per unit length. 

3. The Numerical Scheme 

The numerical scheme consists of a two step, MacCormack-type predictor-corrector methodol- 
ogy [Peyret and Taylor, 1983] for each of the equations. For convenience, the 2D mass, momentum, 
and energy conservation equations can be written in the form 



T = 





(7) 



d t f + ft (F[f]) + 2 (G[f]) = 0, 




where the square parentheses indicate functional dependence, specifically 



F[f\ = F{i,d 1 i,d 2 i,Vp). 



At any grid point the predictor step 




(9) 
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is followed by the corrector step 



2 Axi 
2 As 2 
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where Aj W( i stands for the forward difference, A^ for the backward difference in the relevant 



direction, and V 2 ctr for the center-differenced Laplacian operator. 



4. Numerical Simulations 

There are three dimensionless parameters of interest in the problem: the Reynolds number, 
the Weber number, and a Mach number: 

A= *™ w,= ™ Mc=V (11) 

The number Ma c appears immediately in the non-dimensional form of the pressure, (5). However, 
it is physically more appealing to consider the numbers Ma(p) = V/c s where c 2 s = dpo/dp is the 
sound speed squared, and depends on the density p. 

Towards validating the numerical model, we first present the results of a few test cases. In Fig. 1, 
at a given temperature, the solid line gives the densities of the two coexisting phases calculated 
using the Maxwell equal area construction. The symbols are the densities of the liquid and gas 
phases obtained from the numerical simulations. The simulation consisted of equilibration of a flat 
interface between the high and low density phases on a 64 X 64 grid. The Reynolds number in (II) 
(based on the grid spacing As) was set at 2.0 and the Weber number (also based on As) was set at 
1.0. Using formula (7), the value of the nondimensional surface tension coefficient (ctdi m Ax / (p 2 c X), 
where a. dim is the dimensional surface tension coefficient) is calculated for future reference to be 
a(0/0 c = 0.85, We = 1.0) = 0.47 and a(0/0 c = 0.85, We = 0.6) = 0.64. 



Next we verify Laplace's law for surface tension. The simulations here consisted of letting go to 
equilibrium drops of the liquid phase suspended in the gaseous phase. The computational domain in 
this case was either 256x256 or 64x64 and the Reynolds number was set at 2.0 (again based on As). 
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Fig. 1 The densities of the coexisting phases at different temperatures. The solid line results from 
a thermodynamic calculation in which the Gibbs free energies of the two phases at the given 
temperature are equated. The symbols are results from the numerical simulations. 

The Weber number was either 0.6 or 1.0 and the surface tension coefficient a calculated from the 
flat interface simulations were used to plot a/R on the s-axis and the measured pressure difference 
between the bulk liquid and gas phases Ap is plotted on the y-axis in Fig. 2. The deviation from the 
least squares fit is within 5%. The linearity of the plot verifies Faplace's law and the consistency 
in the model of the surface tension, noting that the surface tension coefficient was obtained from 
the flat interface simulations. We next present two simple applications of the model — a case of 
instability of a high Reynolds number liquid jet in the gaseous phase and an instance of isothermal 
spinodal decomposition. 




Fig. 2 Verification of Laplace's law. On the s-axis we plot ^, where R is the radius of the drop and 
a is the value of the surface tension coefficient calculated from the above simulation of flat 
interfaces. On the y-axis is the numerically measured pressure difference between the inside 
and the outside of the drop. The solid line is the least-squares linear fit and has a slope of 
f .037 with a standard deviation of 0.0008. 

4.1 Instability of a planar jet 

We consider the evolution of a high Reynolds number planar jet of the liquid phase in the 
gaseous phase. We note that there are very few simulations of two-phase flows at high Reynolds 
number considering their computational difficulty. Fig. 3 shows the shape of the jet at times 0, 200, 
400, 800, and fOOO nondimensional units. Initially (t=0), the Reynolds number of the jet based on 
the diameter is about 800, the Weber number, also based on the diameter is about 80, and Ma c is 
about 0.23. 512 grid points are used along the axis of the jet and 128 grid points transverse to it. 
The boundary conditions are doubly periodic and the computations are isothermal. The instability 
of the jet is seen as the unsteady deformation of its shape. Note that in the absence of forcing, 
the energy of the jet is constantly decaying with time. In a second experiment, the same initial 



Fig. 3 Instability of a planar liquid jet in the gas phase at an initial Reynolds number (based on 
the diameter) of about 800 and a Weber number (based on the diameter) of about 80. The 
snapshots of the liquid jet at time 0, 200, 400, 800, and 1000 are shown from top to bottom 
for this isothermal simulation. 

jet configuration is used, but the surface tension is increased by decreasing the Weber number to 
about 40. The snapshots of the jet at the same times as in Fig. 3 are shown for the case with the 
increased surface tension in Fig. 4. The simulation clearly show how the increased surface tension 
tends to stabilize the jet. 



Fig. 4 The same setup as in Fig. 3, but the surface tension is increased by a factor of 2 from that 
case. The snapshots are at the same times as before, and the increased stability of the jet is 
clear. 

4.2 Spinodal Decomposition 

When a solid or fluid is rapidly quenched from above the critical point where there is only 
one homogeneous phase to below the critical point where two or more different phases may coexist 
(within the spinodal), the homogenous phase is unstable and therefore separates spontaneously into 
domains occupied by the coexisting phases. After that, the average size of the domains tends to 
increase in an effort to decrease the interfacial energy. Such spontaneous formation of domains and 
their subsequent coarsening is called spinodal decomposition. The rate of growth of the average size 
of these domains is both of theoretical interest and of practical value. When the system undergoing 
spinodal decomposition is a fluid system, the advection of fluid particles leads to an enhanced growth 
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rate as compared to — say — alloys where there are no advective processes and the growth is limited 
by diffusion. While the precise mechanisms and the growth exponents are not well understood in 
the presence of fluid dynamics, particularly when aspects like compressibility are concerned, it is 
expected that there exist a few universality classes characterizing the growth of domains in such 
systems. This expectation is based in part on the success of scaling theories in conjunction with 
numerical and experimental studies in binary alloys [Gunton, San Miguel, and Sahni, 1983; Rogers, 
Elder, and Desai, 1988; and references therein] and some such similar successes in immiscible binary 
fluid systems [Ma, Maritan, Banavar, & Koplik, 1992; Farrell and Vails, 1990; Furukawa, 1985; and 
references therein]. The numerical study of spinodal decomposition in liquid- vapor systems however 
has not received much attention (in contrast to binary alloy and immiscible binary fluid systems) 
mainly because of a lack of models including all of the fluid dynamics and the energetics (exceptions 
are Langevin fluid models [Farrell and Vails, 1990]). The present compressible model is simple and 
complete except for the effects of temperature variations and the release of latent heat. 

The initial condition is a uniform state at a supercritical temperature with a noise of ampli- 
tude 0.2p c in the density field superposed over an average density of 1.06344p c . The system is then 
quenched to 0.85# c . After this temperature quench, the homogenous phase finds itself in a region of 
mechanical instability (pressure decreases with increasing density) and therefore separates. Fig. 5 
shows the snapshots at times 25, 50, 250 and 500 nondimensional time units of the spinodal decom- 
position following the temperature quench. A doubly periodic 512 X 512 domain was used. At the 
final time shown (t=500) in Fig. 5, Ma c is 1.52, based on the velocity averaged over the domain. 
The Reynolds number in the dense phase, based on the average domain size, is about 95, and the 
Weber number (again based on the average domain size) is about 120. 

By considering the circularly- averaged two-point correlation function (the second-order struc- 
ture function) of the density at a given time, the average domain size is estimated as the first zero 
crossing of the correlation function. The domain growth with time is shown as symbols on a log-log 
plot in Fig. 6. Here the data is obtained by an ensemble average over 12 different runs in which only 
the seed of the random number generator was varied to initialize the noise in the density field. For 
the domain sizes lying between 20Aa; and 128Aa; (the computational domain was 512Aa; X 512Aa;), 
a least-squares linear fit gave us a slope of 0.70 with a standard deviation of 0.01. The data for 
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Fig. 5 Spinodal decomposition in the van der Waals gas following a temperature quench at time 
0. The four snapshots of the density field are at times 25, 50, 250, and 500 nondimensional 
time units. The average density is 1.063 and the density of the coexisting liquid and gas 
phases are respectively 0.319 and 1.804 (across a flat interface). 



which the power law growth was estimated (the late stage) is plotted as ' + ' and the rest of the data 
as open squares. The region of Fig. 6 over which the data was fit with a power law is replotted to 
show the fit more clearly. We plan to present a more complete study of this problem elsewhere, but 
we point out that in this two-dimensional isothermal simulation where inertial effects are important 
(late stage) the 0.70 growth exponent is close to the higher values of the exponents that have been 
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Fig. 6 The domain size as a function of time on a log-log scale. The data is obtained by an ensemble 
average over 12 different runs, one of which is shown in Fig. 5. The runs vary only in the 
seed used for the random number genrator used to initialize the noise in the density field at 
time 0. The solid line is the least-squares linear fit obtained for the data with domain sizes 
lying between 20Ax and 128Aa;. It has a slope of 0.70 with a standard deviation of 0.01. 
The data used to obtain the fit is shown using ' + ' and the rest using open squares. For the 
region of the data fit, in (b), we have plotted D(t)t~ 2 / 3 in the ' + ' symbols and D(t)t~°' 7 in 
the '*' symbols, where D(t) is the average domain size at time t. 



observed [Farrell and Vails, 1990]. 
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Appendix 

We recall the form of the mass, momentum and energy conservation equations for a fluid whose 
free energy is given by the Van der Waals- Cahn-Hilliard model. The mass and momentum equations 
are identical to those in Section 1. In addition there is the energy equation 

d t {pe t ) + dj(pe t Uj) = djiaijUj) - djqj (12) 

where <7 8 j = + + po^ij, e t = e + u 2 /2 is the total specific energy and qj is an unknown 
energy flux. Next, we postulate an extensive (Helmholtz) free energy F in the form 

F = F + ±\(p,T)\Vp\ 2 , (13) 

where Fo is the bulk free energy. We have generalized the approach of the text to allow for a 
dependence of A on density and temperature. ^From thermodynamics, considering p and T as the 
two independent variables of state, the entropy is given by S = —dF/dF, the derivative being 
taken at constant p. The extensive internal energy is pe = E, and is expressed in terms of the free 
energy by 

E = F + FS. 

With these definitions, it may be argued [Dunn and Serrin 1965; Felderhof 1970; Casal and Gouin 
1975; Evans, 1979, Falk 1992] that the energy flux and the stresses must be of the form 

q t = -kd t F + p\V ■ ud iP , (14) 



aff = -Xdtpdjp + - [\ - p—j \Vp\ 2 + P V ■ (XV p) 

(v) 

where A; is a positive transport coefficient, the a\-' are viscous stresses and the buik pressure po is 
given by po = —Fo + pdFo/dp. The argument is in fact a proof in the perfect fluid, no dissipation 
case and rests on the usuai assumption for entropy production in the dissipative case. We may 
chose the free energy to be in the van der Waals form 

F = -pR g T log (j-bj ~ap 2 

so that the equation of state is in the familiar form 

{p+^)(v-b) = R g T. (15) 

The equation of state and the expression of the stresses are now identical to those in the text if we 
fix 9 = F and make A a constant, then rewrite equation (15) in the form (5). 



